#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Oct 22 20:56:19 2023

@author: liqingsimac
"""
import numpy as np
import statsmodels.api as sm

a = np.loadtxt('data10_7_1.txt')
x = a[:,0]
p = a[:,2]/a[:,1]
X = sm.add_constant(x)
y=np.log(p/(1-p))
md = sm.OLS(y, X).fit()
print(md.summary())
b = md.params

import matplotlib.pyplot as plt
fig=plt.figure()
ax=fig.add_subplot(111)
ax.plot(x,p,'bo')
ax.set_xlabel('x')
ax.set_ylabel('p')

xx = np.linspace(start=-20,stop=30,num=51,endpoint=True)
pp = 1/(1+np.exp(-b[0]-b[1]*xx))
ax.plot(xx,pp,'r-')
